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I.  Introduction 


The  formulas  for  the  determination  of  the  coefficients  of  the  spherical 
harmonics  of  the  earth’s  gravitational  potential  require  the  free-air  gravity 
anomalies  to  be  given  on  a  spherei  at  least  at  a  simple  surface,  e.g.,  the 
ellipsoid.  Thus  we  have  to  continue  the  free-air  gravity  anomalies  downward 
to  a  sphere  or  an  ellipsoid. 

The  validity  of  the  analytically  downward  continuation  of  the  free-air 
gravity  anomalies  inside  the  earth  is  guaranteed  by  Range’s  theorem  (Moritz, 
1980,  p.  67).  Of  course  these  gravity  anomalies  are  not  the  original  gravity 
anomalies  inside  the  earth.  The  downward  continuation  gives  a  fictitious 
gravity  anomaly  on  the  ellipsoid  that  generates  a  disturbing  potential  on  and 
outside  the  earth  that  coincides  with  the  original  disturbing  potential  T  on 
and  outside  the  earth. 

Moritz  (1980)  suggested  that  the  free-air  anomalies  be  continued  to  the 
point  level.  We  can  analytically  take  also  the  ellipsoid  as  the  reference 
surface  and  use  this  method  to  continue  the  gravity  anomalies  down  to  the 
ellipsoid. 

Bjerhammar  (1964)  suggested  using  the  Poisson’s  integral  to  continue  the 
free-air  anomalies  downward  to  a  sphere  embedded  inside  the  earth.  This 
problem  was  solved  by  using  the  discrete  technique  and  matrix  formulas. 

Pellinen  (1966)  studied  the  methods  for  the  determination  of  the 
coefficients  of  the  spherical  harmonics  of  the  earth’s  gravitational  potential 
and  he  added  a  term,  which  can  be  easily  transformed  into  terrain  correction, 
to  the  free-air  gravity  anomalies. 

In  this  work  we  carry  out  some  numerical  investigation  with  the  above 
mentioned  methods.  The  terrain  correction  and  the  Molodensky’s  correction 
term  g,  are  computed  on  a  global  basis. 


2.  Analytical  Downward  Continuation  by  Using  Taylor  Series 


We  continue  the  free-air  gravity  anomalies  downward  to  a  point  level  by 
using  the  Taylor  series  given  by  (Moritz,  1980,  p.  378): 


Ag  =  Ag' 


+ 


az* 


(1) 


where  z  =  h-hp  is  the  elevation  difference  with  respect  to  the  computation 
point  P,  Ag'  is  the  gravity  anomaly  on  the  point  level  and  Ag  is  the  gravity 
anomaly  on  the  earth’s  surface. 

By  using  the  up  and  downward  continuation  operators  the  correction 
terms  of  the  gravity  anomalies  at  the  point  level-  can  be  expressed  as  follows 

g,  =  -2L(Ag) 

ga  =  -zL(g,)  -  %z»LL(Ag)  (2) 


with 
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««  =  £  II  ^ 

<s 

where  io  -  2Rsin'V'/2,  i>  is  the  distance  between  the  current  point  and  the 
computation  point  P|  a  is  the  unit  spherei  and  R  is  the  mean  radius  of  the 
earth. 

If  the  ellipsoid  is  chosen  as  the  reference  surface,  in  other  words,  we 
want  to  have  the  reduced  gravity  anomalies  not  on  the  point  level  but  on  the 
ellipsoid,  we  have: 

=  -hL(Ag) 

g,  =  -hL(g,)  -  %h»LL(Ag)  ,  (4) 


where  h  is  the  height  of  the  topography  above  the  ellipsoid. 

Based  on  Runge's  theorem  the  reference  surface  can  be  chosen  arbitrarily  for 
the  downward  continuation.  One  can  show  that  the  downward  continuation  of 
the  gravity  anomalies  to  the  point  level  or  to  the  ellipsoid  is  equivalent  since 
they  generate  the  same  disturbing  potential  T  on  and  outside  the  earth  (cf. 
Sideris,  1987). 

As  a  first  approximation  we  have  the  gravity  anomalies  on  the  ellipsoid: 

Ag*  =  Ag  +  gi  =  Ag  -  hL(Ag)  (5) 

This  is  called  "gradient  solution"  (Moritz,  1966). 

The  relationship  between  the  first  correction  term  of  the  Molodensky’s 
solution  Gt  and  g|  is: 


s.  =  «■  -  f;  II  d.  . 


(6) 


The  first  correction  term  of  the  Molodensky's  solution  is  given  by 


“■  =  I;  II  • 

a 

After  some  reordering,  equation  (6)  can  be  changed  into 

<r 

where  G'  is  the  term  which  was  suggested  by  Pellinen  (1966,  p.  70): 


(7) 


(8) 


S'  =  !;  II  . 

a 

If  we  assume  there  is  a  linear  relationship  between  the  gravity  anomaly 
and  the  elevation: 
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Ag  =  a  +  bh,  b  =  2nV.p  , 


(10) 


where  a  is  a  constant,  k  the  gravitational  constant  and  p  the  density  of  the 
crust,  then  eqn.  (9)  becomes 

G'  =  C  =  I  kpH^  IJ  da  ,  (11) 

a  ° 

where  C  denotes  the  terrain  correction. 

Equation  (8)  can  be  written  in  the  form 

g.  =  C  -  i  kpR"  JJ  da  .  (12) 

a  ° 

Equation  (12)  gives  the  relationship  between  the  first  correction  term  of 
Moritz’s  solution  with  the  terrain  correction.  Moritz  (1966,  p.  106)  shows  that 
the  terrain  correction  plays  a  role  in  the  estimation  of  the  low  degrees  of  the 
spherical  harmonica  of  the  earth’s  gravitational  potential,  using  surface 
gravity  data.  But  the  last  integral  in  (12)  cannot  be  neglected  if  the  high 
degrees  of  the  spherical  harmonics  are  determined. 

Prom  the  equations  (11)  and  (12)  we  get 

gj  =  -kpR»hp  jj  ^  da  .  (12) 

<r 

We  can  get  equation  (13)  directly  from: 

g,  =  -hpL(Ag)  =  -hp  da  (14) 

a 

under  the  assumption  (10). 

The  formulas  (11)  and  (13)  have  the  advantage  that  they  can  be  evaluated 
from  topograghy  only,  no  gravity  anomalies  being  needed.  The  weakness  is 
that  one  has  to  make  the  assumption  (10). 

In  a  plane  approximation,  eqns.  (11)  and  (13)  become 


g,  =  -hp  kp  IJ  dxdy  ,  (15) 

C  =  jkp  JJ  dxdy  (16) 

T 

X 

with  t  =  [(x-x  )*  +  (y-y  )*]^  •  where  r  is  the  two-dimensional  plane. 
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3.  Downward  Continuation  of  the  Free-Air  Gravity  Anomalies  Using  the 
Poisson's  Integral 

Another  method  for  the  downward  continuation  of  the  gravity  anomalies  to 
the  ellipsoid  is  using  the  Poisson’s  Integral  (Moritz,  1966,  sec.  7).  In  the 
plane  approximations  we  have: 

r 

where  Ag*,  Ag  are  the  gravity  anomalies  on  the  ellipsoid  and  on  the  earth’s 
surface  respectively.  The  distance  between  the  computation  point  P  on  the 


earth’s  surface  and  the  current  point  on  the  ellipsoid  is 

t  -  [(x-Xp)*  +  (y-yp)*  +  hp*]^  .  (18) 

If  we  introduce  the  identical  operator  I  and  the  operator  Lp 

If  =  f  (19) 

V  =  k  11  ^  •  (20) 

T 

Then  eqn.  (17)  can  be  written  in  the  form 

(I  +  Lp)Ag»  =  Ag  .  (21) 

Its  solution  is  given  by 
Ag*  =  (I  +  Lp)"»Ag 

=  7^(-l)"Lp"  Ag  (22) 


n'O 

<D 

0 

with 

g*  =  (-l)"L^Ag  .  (23) 

More  explicitly  we  write 
g*  =  Ag 

gf  =  -LpAg  =  -  11  dxdy 

T 

gf  -  Lp*Ag  =  -  Lpg^c  =  -  11  dxdy  (24) 


Comparing  gf  in  (24)  with  (14),  we  find  that  both  formulas  are  similar.  The 
difference  is  that  the  integral  kernel  for  g|  is  and  for  g,  is  <“’.  It  was 

shown  (Moritz,  1966,  p.  60)  that  the  formulas  for  gf  and  g,  are  approximately 
equivalent,  gf  will  give  more  accurate  results. 
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L  Some  Considerations  About  Applying  the  FFT  Technique  to  the  Poisson's 
Integral 

The  PFT  (Fast  Fourier  Transformation)  is  widely  used  for  numerical 
computations.  We  consider  here  how  the  FFT  technique  can  be  used  for 
P-integral  (Poisson’s  integral).  This  work  has  been  described  by  Euler  et  al. 
(1980),  but  we  have  to  recognize  that  the  Fourier  transformation  can  not  be 
directly  used  for  the  P-integral  in  plane  approximation  due  to  the  elevation 
variable  which  is  included  in  the  kernel  of  the  P-integral.  We  have  to  do 
some  things  before  we  apply  the  FFT  to  the  P-integral. 

In  plane  approximation  the  P-integral  becomes 

=  Is  fl  7?  *“>>’  <===' 

T 

where  Agp  is  the  gravity  anomaly  at  point  P  which  is  located  on  the  earth’s 
surface,  4g*  is  the  gravity  anomaly  on  the  ellipsoid,  hp  is  the  elevation  of  the 
point  P.  The  distance  between  point  P  and  the  current  point  on  the  ellipsoid 
is 

t  =  '/U-Xp)^+(y-yp^^+hJ  (26) 

Because  the  elevation  variable,  hp,  is  included  in  the  kernel  i~^ ,  the  kernel 
can  not  be  written  in  the  form 


<  =  <(Xp,yp,x,y,h(Xp,yp))  =  ^(x-Xp,y-yp) 

and  the  P-integral  is  not  a  convolution.  The  definition  of  the  two  dimensional 
convolution  is  given  in  Bracewell  (1965,  p.  243). 

The  two  dimensional  Fourier  transformation  and  its  inverse  are  defined  as: 


F(u,v)  =  F(f(x,y)}  = 

f(x,y)  =  >-‘{F(u,v)} 
where  the  symbols  F,  F~ 


iJ  f(x,y)e  (xu+avljjxdy 
=  iJ  F(u,v)e^”^’ (’‘“■’■y^'^dudv 


(27) 


(28) 


denote  the  Fourier  transformation  and  its  inverse. 


If  the  elevation  hp  is  constant  in  the  integration  region,  the  integral  (25) 
can  be  considered  as  a  convolution.  Apply  the  Fourier  transformation  to  (25) 
and  we  get  immediately  (Euler,  et.  al.,  1986) 

G(u,v)  =  e“^^^’P“G*(u,v)  (29) 

with  ui  =  Vu^+v*  ,  u,v  are  the  frequency  variables,  G  and  G»  are  the  Fourier 
transformations  of  the  gravity  anomalies  Ag,  4g*  respectively. 


From  (29)  we  have 

G*(u,v)  =  e^"^P“G(u,v) 
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(30) 


(30)  shows  that  the  gravity  anomaly  on  the  ellipsoid  is  rougher  than  the 
gravity  anomaly  on  the  earth’s  surface.  Although  the  relationship  of  the 
spectra  of  the  gravity  anomalies  dg  and  dg*  is  really  not  so  simple  as  (30) 
shows,  the  relationship  of  the  spectra  of  the  gravity  anomalies  should  have 
similar  properties;  the  upward  continuation  makes  .  the  gravity  anomalies 
smoother  and  the  downward  continuation  makes  the  gravity  anomalies  rougher. 
We  have  to  see  that  if  we  have  some  observation  errors  in  the  gravity 
anomalies,  and  normally  the  errors  are  assumed  to  be  composed  of  high 
frequencies,  e.g.,  the  white  noise,  then  the  downward  continuation  enlarges 
the  spectra  of  the  errors  by  the  factor  exp(27rhpCi>).  This  means,  a  small  error 
in  dg  can  cause  a  big  error  in  dg*,  if  the  errors  in  dg  are  composed  of  high 
frequencies.  Therefore  the  downward  continuation  may  not  be  stable.  The 
situation  is  not  so  serious  in  practice  if  we  use  the  mean  block  values  which 
may  decrease  the  effect  of  the  high  frequencies  in  the  errors. 

Now  we  go  back  to  see  how  we  can  apply  the  Fourier  transformation  to 
the  P-integral. 

For  the  Kernel  function  we  have 

=  (x-xp)»  +  (y-yp)*  +  hp» 

=  (x-xp)*  +(y-yp)*  +  C  +  (hp»  -  c») 

=  S"  +  (hp*  -  C*)  ,  (31) 

with 

S*  =  (x-Xp)»  +  (y-yp)*  +  C*  , 

C  is  a  constant.  If  we  take  C  =  hn,ax>  maximum  of  the  elevation  in  the 

integration  region,  we  have 

t2  =  s*(l  +  )  =  S^(l  +  fl  )  .  (32) 


with 


H  = 


hp*  - 


h^ax 


The  P-integral  is  expanded  into  the  series 


dgp  = 


h^i 

27T 


f  Aa* 

)} 


3  H 


15  H’  105 


/  ^  It _  X  IX  \  1  J 

9  Q  04  AO  dxdy 


Because  we  have  always 


the  series  in  the  integral  (33)  converges  absolutely. 
Now  we  apply  the  Fourier  transformation  to  (33) 


(33) 


(34) 
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For  the  kernel  S  we  have 

2  2  ,  .  2 


S-(2n+0  .  (-  £  )(_  I 


an 


2n-l  ^  a(h^„)n 


(  i  ) 


n=l , 2, . . . 

The  Fourier  transform  of  the  kernel  S“*  is  (Dracewell,  1965,  p.249): 

J  i  g-2Trhfflaxul 

S  u< 


(35) 


(36) 


By  using  eqns.  (35)  and  (36)  we  have 


^  s-  ^  2  a(h^„J 


S»  ^  3  3(h-,3J 


£T  .1.  'V  =  A~3nhmaxui 


.  nl}= 


max  ^max 


(r— ^  +  27tuj) 


(37) 

(38) 


then  (34)  becomes 

Ago  =  ho[/^'{G*  T-^  —  F^*{G*(ir^  +  2nui)  +  ...] 

"  ^ax  “max 

(39) 


Equations  (34)  and  (39)  have  more  theoretical  meaning  because  the  series 
converges  slowly.  Numerical  tests  show  that  if  a  1  mgal  accuracy  of  the 
P-integral  is  needed,  the  computation  should  include  up  to  the  fifth  term  in 
the  series  (34). 


In  order  to  increasse  the  convergence  of  the  series  we  modify  this 


method. 

We  decompose  the  integral  kernel  t  ’  into  two  parts: 

=  £73  + 

d— 3 

(40) 

with 

=  ( 

1  0 

x,y  cD 
x,y 

(41) 

=  (  ° 

1  £-3 

x,y  cD 
x.y  ^ 

(42) 

where  D 

is  the 

innermost  zone  surrounding  the  point  P. 

Now  the  P-integral  is  decomposed  into  two  integrals: 


Agt 


h^ 

2n 


if 


II  if  *“*>' 


r  T 

We  expand  the  second  integral  in  the  series  as 


ha 

2n 


ffr  =  I: 


If 


Atf*  ,,  3  H  ^  15  . 

S!i  ^  2  S?i  8  Sh 


(43) 


(44) 
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with 


^ii 


x,y  i»D 
x,y  sD 


(45) 


The  integrals  in  (44)  are  not  difficult  to  compute  by  the  use  of  the  Fourier 
transformation  or  FFT  technique.  If  the  innermost  zone  D  is  properly  chosen, 
the  series  (44)  converges  very  fast.  We  give  a  numerical  teat  hereafter. 


For  the  first  integral  in  (43)  the  numerical  integration  is  used.  Because 
the  integral  is  only  evaluated  in  the  innermost  zone,  the  computation  can  be 
completed  fast. 

Now  we  give  a  numerical  test  for  the  above  mentioned  methods. 

We  simulate  the  gravity  anomalies  data  on  the  ellipsoid  by  using  the 
formula 


Ag*  =  O.llh  -  200.  (46) 

where  h  is  in  aeter;  ig*  in  mgals. 

The  digital  elevation  at  30"  grid  interval  for  the  USA  is  used.  The  extent 
of  the  test  region  is  50  x  50  km.  In  this  region  the  maximum  elevation  is  3260 
meters;  the  lowest  elevation  is  1830  meters  and  the  mean  elevation  is  2374 
meters.  Values  have  been  computed  for  points  located  in  a  rectangle  in  the 
center  of  the  area  where  the  points  are  approximately  I  km  apart. 

Results  for  Agp  from  (34)  up  to  the  fourth  term  are  given  in  Table  1.  The 
numerical  integration  values  of  the  P-integral  are  also  given  as  the  check 
values. 


Table  1.  Results  for  Agp  from  (34)  (Unit:  mgal) 


Point  number 

1 

2 

3 

4 

1st  term 

41.11 

41.64 

38.71 

37.89 

2nd  term 

10.48 

10.43 

11.01 

11.20 

3rd  term 

3.53 

3.44 

4.16 

4.39 

4th  terra 

1.32 

1.28 

1.74 

1.90 

Sum  (Agp) 

56.44 

56.77 

55.62 

55.38 

Check  Values 

57.25 

57.53 

56.68 

56.78 

Sum-Check  Val. 

-0.81 

-0.76 

-1.06 

-1.40 

Ag* 

66.20 

68.40 

53.00 

47.50 

Table  1  shows  the  series  (34)  numerically  converges  with  factor  1/3 
approximately  in  the  test  region.  If  we  want  to  have  the  computing  accuracy 
at  1  mgal,  we  have  to  complete  the  computations  up  to  the  5th  term  of  the 
series. 

We  now  give  the  results  for  Agp  from  (43)  and  (44)  in  Table  2.  The 
innermost  zone  is  chosen  as  a  5’x5’  square  and  the  P-integral  is  completed  in 
this  square  with  numerical  integration.  The  FFT  technique  is  applied  to  (44). 
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Table  2.  Results  for  Agp  from  (43)  and  (44).  (Unitimgal) 


Point  number 

1 

2 

3 

4 

Innermost  Zone 

39.91 

40.05 

40.02 

40.40 

1st  term  (eq.(44)) 

16.25 

16.38 

15.71 

15.32 

2nd  term 

1.07 

1.06 

1.17 

1.19 

3rd  term 

0.08 

0.08 

0.10 

0.11 

4th  term 

0.01 

0.01 

0.01 

0.01 

Sum 

57.32 

57.58 

57.01 

57.03 

Check  Values 

57.25 

57.53 

56.68 

56.68 

Sum-check  Values 

0.07 

0.05 

0.33 

0.25 

Ag* 

66.20 

68.40 

53.00 

47.50 

Table  2  shows  the  series  (44)  converges  very  fast.  Even  when  we  choose 
a  small  innermost  zone  (5*x5*  square),  the  computations  for  the  series  (44) 
need  only  the  first  and  second  term. 

We  give  a  short  summary  about  the  use  of  FFT  for  the  P-integral.  We  can 
use  FFT  through  the  computation  of  (34).  The  disadvantage  is  that  the  series 
converges  slowly.  The  numerical  teat  shows  if  we  want  to  have  computing 
accuracy  at  1  mgal,  we  have  to  compute  the  series  (34)  up  to  the  5th  term. 
In  the  modified  method  we  separate  the  integration  region  into  two  parts: 
innermost  zone  and  outer  zone.  We  evaluate  the  P-integral  in  the  innermost 
zone  with  numerical  integration.  The  innermost  zone  is  chosen  small  and  this 
computation  is  not  time  consuming.  Equation  (44)  is  used  for  computing  the 
P-integral  in  the  outer  zone.  For  the  numerical  computation,  eqns.  (43)  and 
(44)  are  more  beneficial  than  (34). 

In  these  procedures  there  are  several  concerns.  First,  we  must 
understand  why  FFT  calculations  near  the  boundaries  of  the  data  area  give 
incorrect  results  as  has  been  demonstrated  by  Sideris  (1985,  p.  70,  Figure  1). 
In  addition  we  must  consider  how  mean  block  values  can  be  used  in  the  FFT 
computations. 

The  integrals  in  physical  geodesy  are  defined  on  a  sphere  or  approximated 
in  a  plane.  In  practice  the  integrals  are  evaluated  in  a  limited  region. 
Generally  we  write 

a  b 

g(xp,yp)  =11  f(x,y)  k(x,y,Xp,yp)  dxdy  (47) 

O  O 

where  k(x,y,Xp,yp)  is  the  integral  kernel.  If  the  kernel  satisfies  the 
condition 

k(x,y,Xp,yp)  =  k(x-Xp,y-yp)  (48) 
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r 


then  the  discrete  Fourier  transformation  can  be  used  for  (47).  Here  we 
consider  only  this  case. 

In  order  to  use  the  FFT  technique,  we  first  have  to  discretize  the  integral 
in  (47)  in  the  form: 


NI-1  N2-2 

gp  =  ^  ^  '  f(mAx,nAy)  k(inAx-Xp,nAy-yp)  AxAy  (49) 

in=  o  n-  o 


where  Ax, Ay  are  the  variable  increments,  Nl,  N2  are  the  grid  point  numbers  in 
the  integration  region  along  x,y  directions. 

Setting 

(  Xp  =  rAx  r  =  0,1,...N1-1 

I  yp  =  sAy  s  =  0,1,-..N2-1 

(49)  becomes 

Nl-1  N2-1 

gp  =  J  '  f(mAx,nAy)  k[  (m-r)Ax,  (n-s)  AyJ  AxAy  (51) 

m=o  n=o 


If  we  want  to  apply  FFT  to  (51),  we  have  to  extend  the  function  f  and  the 
kernel  k  into  the  whole  plane  as  follows: 


Figure  1.  Periodic  extension  of  the  function  f  and  kernel  k, 
This  periodic  extension  means  we  now  have 
f [ (m+Nl) Ax, (n+N2) Ay)  =  f(mAx,  nAy) 
k)  fnrt-Nl) Ax,  (n+N2) Ay]  =  k(mAx,nAy) 


(y  -  0) 


(52) 
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Here  we  must  recognize  that  this  extension  causes  errors  in  the  computations 
because  the  function  f  and  the  kernel  k  are  not  periodic. 


The  discrete  Fourier  transformation  and  its  inverse  are  defined  as 


F(mdu,nAv) 


Ml-1  N2-1  .  . 

'  f(inAx,nAy)e*"*  H2^ 

k-o  t=o 


(53) 


f (kAx, lAy) 


_ 1_ 

N1N2 


Nl-1  N2-1 

nz  Ez 


n»=o  n-o 


F(inAu,nAv)e  m?)  . 


(54) 


The  discrete  Fourier  transformation  of  (51)  is 

G(mAu,nAv)  =  ^  '  gp(kAx, iAy)e*^‘ (^ 

k.A 

^  ^  f(rAx,sAy)k((r-k)Ax,  (s-A)Ay]e*"‘ (^  "^  (55) 

k,i  r,s 


Setting 


fr-k=g  g  =  -Nl+A+r, . . . ,r 

ls-i=h  h=  -N2+<+s, . . . ,3 


(56) 


(55)  becomes 


G(mAu,nAv)=  ^  '  f  (rAx,sAy)e*’^‘ (™ 

r.s 

•  ^  ^  k(gAx,hAy)e”*^‘ ^)AxAy  (57) 

k,i 

=  F(mAu,nAv)K(mAu,nAv)AxAy  , 


where 

F(mAu,nAv)=  ^  ^  '  f  (rAx,sAy)e*"’ (™  ■*"  »  (^^^ 

r=o  s=o 

K(mAu,nAv)=  ^  |  k(gAx,hAy)e“’’^‘ 

k.i 


r-Nl+1  S-N2+1 


-nzi 


k(gAx,hAy)e 


(59) 


«=r 


h=3 


If  the  kernel  is  periodic  as  (52)  shows,  we  have  immediately 
Nl-1  N2-1  , 

K(niAu,nAv)  =  J  '  j  '  k(gAx.hAy)e*"‘ .  (60) 

g=o  h=o 

It  is  not  the  case  in  the  numerical  integration  because  the  kernel  function  is 
not  periodic.  In  order  to  see  this  more  specifically,  we  give  an  example.  We 

take  the  kernel 

k(x,y)  =  s“’ 

=  (x*  +  y»  +  h^a*)*/* 

Figure  2  represents  the  kernel  k  and  its  periodic  extension. 


Figure  2.  Kernel  k  and  its  periodic  extension 

The  striped  part  is  the  range  of  the  kernel  s"’  in  (59).  If  we  use  (60),  then 
the  range  of  the  kernel  8~^  is  the  striped  part  plus  the  shaded  part.  The 
shaded  part  is  the  error  caused  by  making  the  kernel  function  periodic.  At 
r=0,  corresponding  to  the  boundary  point,  the  errors  will  go  to  the  largest 
(see  Figure  2).  Therefore  the  boundary  values  of  the  FFT  computations  are 
wrong. 

We  can  see  that  only  at  the  center  of  the  integration  region,  where 


do  we  have  equation  (59)  equal  to  (60).  The  closer  the  computation  points  are 
to  the  center,  the  smaller  the  shaded  part  in  Figure  2,  or,  the  smaller  the 
error  caused  by  extending  the  kernel  function  periodically. 


12 


In  the  FFT  computations  we  use  (60)  instead  of  (59).  Based  on  the  above 
discussion  we  know  only  at  the  center  point  there  are  no  errors  caused  by 
extending  the  kernel  function  periodic.  The  situation  is  not  so  serious  if  the 
kernel  function  decreases  very  fast,  e.g.,io”’  or  ***  case,  the 

errors  are  small  if  the  computation  points  are  not  at  the  boundary  region  (cf. 
Sideris,  1985,p.  70,  Fig.  1).  What  we  have  to  remember  is  that  we  take  a 
proper  boundary  for  the  computations  in  order  to  guarantee  the  computations 
are  sufficiently  accurate. 

This  problem  seems  more  serious  for  applying  FFT  to  the  Stokes  integral, 
because  the  kernel  io*  decreases  slowly.  If  the  boundary  is  not  correctly 
taken,  the  errors  may  be  significant.  In  other  words,  the  boundary  has  to  be 
chosen  large  enough. 

In  order  to  decrease  the  errors  caused  by  extending  the  kernel  function 
periodic,  we  can  enlarge  the  integration  region  by  putting  a  zero 
boundary — some  zero  values  are  assigned  to  the  function  f  in  the  added 
boundary.  Now  we  consider  how  we  can  use  the  block  mean  values  for  the 
FFT  computations. 

If  the  function  f  is  given  in  mean  block  values,  then  the  integral  (47) 
becomes 


m  a 

g(xp,yp)  =  {  I  f(x,y)k(x-Xp,y-yp)dxdy 

O  O 

=  j  '  f(iAx,jAy)  I  I  k(x-Xp,y-yp)dxdy 

»  11  I  i  k  '  1 


where  f(iAx,jAy)  is  the  mean  block  values  of  the  function  f  in  the  block  Aij. 
We  denote  the  integral  in  (61)  by  c  and  put  (50)  in  this  integral 

c[(i-r)Ax, ( j-s)Ay]  =  |  |  k(x-Xp,y-yp)dxdy 

Aij 

(i+7)Ax  (j-Hj)Ay 


I  k(x-rAx,y-sAy)  dxdy 


then  (61)  becomes 


g(rAx,sAy)  = 


c 


(i-j)Ax  (j-j)Ay 


f(iAx,jAy)  c[(i-r)Ax, (j-s)Ay] 


The  discrete  Fourier  transformation  of  (63)  is 

G(kAu,lAv)  =  F(kAu,lAv)  C(kAu,lAv)  (64) 

where  F,  C  are  the  discrete  Fourier  transforms  of  the  mean  values  f  and  the 
kernel  function  c,  respectively. 
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Taking  the  inverse  of  the  discrete  Fourier  transformation  of  eqn. 
( 64 ) ,  we  get 

ffl-1  W2-I  .  ,  , 

g(r4x,siiy)  =  ^  ^  G(kAu, iAv)e-a^<  +  -if) 


k=o  1=0 


Nl-1  N2-1 

ik 


F(kAu, AAv)C(kAu, iAv)e 


k=o  1=0 


We  point  out  that  we  consider  the  FFT  technique  only  as  a  computing  tool. 
With  its  help  we  can  evaluate  (63)  very  fast.  Equations  (64)  and  (65)  give  the 
same  results  as  (63),  except  for  some  errors  caused  by  extending  the  kernel 
function  periodically. 


Numerical  Investigation  of  the  Downward  Continuation  -  A  Comnarison  of 


Downward  Continuation  by  Using  the  Gradient  Solution  and  Terrain 


Correction 


The  downward  continuation  of  the  gravity  anomalies  to  the  ellipsoid  by 
using  the  gradient  solution,  Poisson's  integral  and  the  terrain  correction 
require  a  dense  distribution  of  gravity  anomaly  data  and  elevation  data 
surrounding  the  computation  point.  If  we  use  the  mean  block  values  in  the 
integration,  the  small  block  size  will  be  needed  for  a  high  computational 
accuracy.  Up  to  date  the  mean  free-air  anomaly  in  30'  block  and  the  mean 
elevation  in  5'  block  are  given  in  a  global  basis.  The  mean  free-air  anomaly 
in  30'  block  is  not  suitable  for  the  above  mentioned  integrals. 

The  digital  elevation  at  30"  grid  interval  for  the  U.S.  is  available.  We  will 
use  this  data  for  the  computation  of  (15)  and  (16).  From  the  results  the  mean 
block  values  of  the  gradient  solution  g|  and  the  terrain  correction  (TC)  in  5’ 
and  30'  block  sizes  will  be  formed.  The  5'x5'  mean  elevation  will  be  formed  by 
averaging  the  30”  point  elevation  and  are  used  for  the  computations  of  (15) 
and  (16).  We  want  to  see  what  happens  if  the  bigger  mean  block  values  are 
used. 

The  tests  are  taken  in  three  different  topographic  regions.  In  the  region 
36*<4«37*,  241‘’<\^242*  the  topography  varies  from  200  meters  to  3600  meters; 
in  the  region  39"<4<40’,  253*<X<254’  the  topography  ranges  from  2100  meters 
to  3100  meters;  in  the  region  35*<4*36',  277'«X^278*  the  topography  varies 
from  250  meters  to  1400  meters. 


The  gi  term  and  the  terrain  correction  are  given  in  Table  3-8  for  the  test 
regions.  The  maximum  of  the  gt  term  in  a  5'x5'  mean  block  is  70  mgals  in 
the  rough  area  (36'<4<37*,  241*<X<242’),  and  the  minimum  of  the  g,  term  is 
-28  mgals  (39*<4<40*,  253*<X^254*).  The  maximum  of  the  terrain  correction  in 
a  5'x5'  mean  block  is  22  mgals  (36*^4^37*,  241*<X<242*). 
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Table  3.  gi  in  5’,  30’  and  1**  Mean  Block  Values 


(30"  Point  Elevation  Data) 


1  10 

15 

m 

10 

30 ' x30 '  Mean  Val ues 


10 


l*xl*  Mean  Value 


8 

8 

11 

17 

15 

11 

11 

19 

11 

7 

10 

13 

19 

18 

19 

13 

10 

14 

11 

16 

14 

7 

7 

12 

8 

9 

13 

12 

14 

17 

13 

16 

13 

7 

6 

14 

7 

7 

7 

7 

8 

11 

11 

14 

19 

8 

5 

11 

6 

9 

7 

11 

11 

11 

11 

9 

20 

13 

6 

6 

6 

10 

12 

14 

14 

11 

12 

11 

14 

22 

12 

6 

9 

11 

13 

12 

10 

9 

10 

6 

9 

16 

5 

11 

12 

11 

11 

8 

4 

5 

13 

11 

4 

12 

10 

6 

10 

9 

6 

•  4 

13 

11 

4 

7 

9 

10 

5 

5 

8 

5 

3 

9 

4 

9 

11 

9 

6 

6 

m 

4 

3 

10 

4 

7 

8 

8 

6 

8 

mi 

4 

3 

9 

a)  5’x5’  Mean  Values 


Table  4.  TC  in  5’,  30’  and  1*  Mean  Block  Values 
(30"  Point  Elevation  Data) 
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b)  30'x30'  Mean  Values 


c)  1‘xl*  Mean  Value 


Table  5.  gt  in  5’,  30’  and  1*  Mean  Block  Values 
(30"  Point  Elevation  Data) 
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b)  30' X  30'  Mean  Values 
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c)  I'x  1*  Mean  Value 


a)  5’x  5*  Mean  Values 


Table  8.  TC  in  5’  ,  30’  and  1*  Mean  Block  Values 
(30"  Point  Elevation  Data) 


From  the  above  Tables,  we  see  that  the  correction  term  gi  has  large 
differences  with  the  terrain  corrections  in  the  5’x  5*  mean  block  values.  The 
results  in  bigger  mean  block  values  agree  betler.  The  mean  block  values  of 
g,  agree  very  well  with  the  mean  terrain  correction  in  l*xl*  block. 

We  also  use  the  mean  elevation  in  5’x5’  block  formed  by  averaging  the  30" 
point  elevation  data  for  the  computations  of  (15)  and  (16).  The  results  are 
given  in  Tables  9  and  10. 
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a)  5’x  5’  Mean  Values 

Table  9.  gi  in  5’,  30’  and  1*  Mean  Block  Values 
(5’  Mean  Elevation  Data) 


b)  30' X  30'  Mean  Values 
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c)  l*x  1*  Mean  Value 
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30  Mean  VaJ  ues 


1  c)  I'x  1"  Mean  Value 
1 

1  106* 


a)  5’x5’  Mean  Values 

Table  10.  TC  in  5*,  30’  and  1*  Mean  Block  Values 
(5’  Mean  Elevation  Data) 


Comparing  Table  9  with  Table  5,  we  find  that  if  we  use  the  5’x5’  mean 
elevation  data,  the  mean  block  values  of  gi  will  be  smoother  and  the 
magnitude  becomes  smaller.  For  the  terrain  correction  we  have  the  similar 
conclusion  (cf.  Table  10  and  Table  6).  Since  the  terrain  correction  depends  on 
the  roughness  of  the  topography,  the  results  of  the  terrain  correction  are 
sensitive  to  the  smoothing  effects  (using  the  mean  block  values  is  a  smoothing 
in  the  data). 

In  the  following,  we  want  to  see  what  the  results  will  be  if  we  interpolate 
the  data  to  small  grid  intervals.  We  will  use  the  5’x5’  mean  elevations  and 
consider  this  data  as  point  values  located  at  the  center  of  the  blocks.  Then 
we  interpolate  the  elevation  on  to  a  I’xl’  grid  by  using  a  bicubic  spline 
function.  The  results  are  given  in  Tables  11,  12,  13  and  14. 
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a)  5’x  5’  Mean  Values 

Table  12.  gj  in  5*x5’,  30*x30'  and  l*xl"  Mean  Block 
Values  (5’x5’  Mean  Elevation  Data; 
no  Interpolation) 
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a)  5’x5’  Mean  Values 

Table  13.  TC  in  5’x5’,  30’x30*  and  I’xl*  Mean  Block 
Values  (5’x5’  Mean  Elevation  Data; 
Interpolated  in  I’xl’  Mean  Elevation) 
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30 'x  30'  Mean  Values 
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I'x  1*  Mean  Value 


a)  5’x5’  Mean  Values 

Table  14.  TC  in  5’x5’,  30’x30’  and  l*xl"  Mean  Block 
Values  (5’x5’  Mean  Elevation  Data; 
no  Interpolation) 

Tables  11-14  show  that  there  6U‘e  no  significant  improvements  in  the 
results  by  using  interpolation  of  the  elevation  to  dense  point  values.  From 
the  computing  process  we  know:  the  integration  in  the  innermost  zone  (5’x5' 
square)  has  no  significant  contributions  to  g|  and  TC,  because  we  set  the 
integration  in  the  innermost  zone  zero  when  there  are  no  interpolations  used. 

The  gradient  correction  terms  for  the  gravity  anomalies  were  also 
computed  by  using  the  potential  coefficient  model  OSU86F  (Rapp  and  Cruz, 
1986).  We  use  the  program  written  by  H.  H.  Rapp  for  the  correction  terms  of 
gravity  anomalies  in  the  test  region  and  compare  the  results  in  Table  15. 

Table  15.  Comparison  of  the  Gradient  Solution,  Terrain  Correction 

and  the  Correction  Terms  From  the  Global  Gravity  Model 
0SU86F  (UniUmgal) 


♦ 

X 

height 

-*4^ 

ih»i^ 

*  ah» 

TC 

253.25 

2632 

-0.04 

253.75 

2707 

0.48 

0.05 

39.75 

253,25 

2721 

-0.62 

-0.08 

39.75 

253.75 

2589 

1.53 

0.04 

1 

35.25 

277.25 

839 

0.27 

35.25 

277.75 

441 

-0.25 

0.00 

35.25 

277,25 

789 

0.68 

0.01 

35.25 

277.75 

848 

0.34 

0.01 

0 

36.25 

241.25 

1406 

1.11 

0.04 

5 

36.25 

241.75 

2400 

4.72 

0.39 

■■ 

36.75 

241.25 

2252 

2.92 

36.75 

241.75 

2352 

3.61 

0.28 

12 

HH 

Table  15  shows  that  the  correction  terms  of  the  gravity  anomalies  by 
using  the  Gravity  Model  OSU86F  are  much  smaller  than  gi  and  TC, 
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6.  Computations  of  the  Terrain  Correction  and  the  Gradient  Solution  on  a 
Global  Baals. 

The  elevation  data  in  5'x  5'  mean  block  values  are  available  from  the 

National  Geophysical  Data  Center,  Boulder,  Colorado.  The  elevation  data  has 

been  used  for  the  computation  of  the  gi  term  and  the  terrain  correction  on  a 
global  basis. 

For  the  computation  of  the  gt  term,  we  make  the  assumption 

&g  =  a  bh  ,  b  =  27rkp 

The  parameters  a  and  b  have  been  taken  as  constants  with  b  equal  to  0.11 
mgal/ meter. 

The  integration  region  is  taken  as  15'  in  latilutde  extent  and  30"  in 

longitude.  The  boundary  (or  overlap)  of  the  integration  is  taken  by  50km. 

This  satisfies  the  accuracy  for  most  situations  (Nod,  1980).  But  the  boundary 
is  not  large  enough  for  the  Himalaya  Mountains.  Figure  3  shows  the 
differences  of  the  gi  term  by  using  a  50km  boundary  and  a  250km  boundary. 
The  maximum  difference  is  12  mgals.  Figure  4  shows  the  differences  of  the 
terrain  correction  by  using  a  50km  and  a  250km  boundary.  The  maximum 
difference  is  6mgals.  In  a  flat  area  the  differences  of  the  gi  term  and  the 
terrain  correction  by  using  the  different  boundaries  are  smaller  than  Imgal 
(see  Fig.  5).  Therefore,  we  have  to  take  larger  boundary  if  the  topography  is 
rough  and  high  in  the  integration  region.  In  the  computations  we  have  taken 
the  250km  boundary  for  the  mountain  areas  and  a  50km  boundary  for  the  flat 
areas. 
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FiKur«s  3.  Differencow  of  tho  Ki  Term  Ueinf^  50  km  and  250  km  HournlarioH  in 
an  Area  Conluiriing  the  Himalaya  Muuniaina; 

(Contour  Interval:  2  mgal) 


LONG  I lUDE 
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Oifferunceu  of  th«  Terrain  Correction  lfnirif(  50  km  an<i  250  k 
Boundaries  in  an  Area  Containing  the  Himalaya  Mountains; 

(Contour  Interval:  1  iiignl) 


LONG! ruDE 


QiTIIibl 


Figure  5.  Differences  of  the  g,  Term  due  to  the  IJ 
Boundaries  in  a  'J’uiMJgraphically  Flat  Area; 
(Contour  Interval:  0.4  mgal) 


The  g,  term  and  the  terrain  correction  are  computed  in  5’  mean  block  values 
by  using  the  FFT  technique.  All  computations  are  completed  at  the 
Instruction  and  Research  Compututer  Center  at  The  Ohio  State  University. 
The  CPU  time  required  was  about  10~’  second  for  each  point  on  the  IBM  3081 
at  OSU. 

Tables  16  and  17  show  the  distributions  of  the  gj  term  in  5’  and  30'  mean 
block  values  according  to  the  magnitude.  gi  are  almost  concentrated  in  the 
interval  [-5,5].  There  are  only  a  few  larger  values  for  the  g,  term. 


Table  16.  The  Distribution  of  the  gi  Term  According  to  the  Magnitude 
(5'  Mean  Block  Values) 


1  Interval  (mgals) 

Numbers 

Percentage 

-85. 

.  -75. 

4 

0.000  % 

-75. 

.  -65. 

24 

0.000 

-65. 

.  -55. 

36 

0.000 

-55. 

.  -45. 

86 

0.001 

-45. 

.  -35. 

274 

0.003 

-35. 

.  -25. 

870 

0.009 

-25. 

.  -15. 

3771 

0.040 

-15. 

.  -  5. 

38846 

0.416 

-  5. 

.  5. 

9123185 

97.771 

5. 

.  15. 

121068 

1.297 

15. 

.  25. 

25428 

0.273 

25. 

.  35. 

7173 

0.077 

35. 

.  45. 

2939 

0.031 

45. 

.  55. 

1391 

0.015 

55. 

.  65. 

5083 

0.054 

65. 

.  75. 

410 

0.004 

75. 

.  85. 

217 

0.002 

85. 

.  95. 

116 

0.001 

95. 

.  105. 

77 

0.000 

105. 

.  115. 

51 

0.000 

115. 

.  125. 

48 

0.000 

125. 

.  135. 

29 

0.000 

135. 

.  145. 

22 

0.000 

145. 

.  155. 

16 

0.000 

155. 

.  165. 

5 

0.000 

165. 

.  175. 

10 

0.000 

175. 

.  185. 

9 

0.000 

185. 

.  195. 

3 

0.000 

195. 

.  205. 

3 

0.000 

205. 

.  215. 

1 

0.000 

215. 

.  225. 

0 

0.000 

225. 

.  235. 

0 

0.000 

235. 

.  245. 

1 

0.000 

25 


Table  17.  The  Distribution  of  the  g,  term  According  to  the  Magnitude 
(30'  Mean  Block  Values) 


Interval  (mgals) 

Numbers 

Percentage 

-15.  .  -  5. 

86 

0.033  X 

-  5.  .  5. 

255893 

98.724 

5.  .  15. 

2198 

0.848 

15.  .  25. 

969 

0.374 

25.  .  35. 

49 

0.018 

35.  .  45. 

4 

0.002 

45.  .  55. 

1 

0.040 

The  terrain  corrections  are  concentrated  between  0  and  10  mgals.  The 
distribution  of  the  terrain  correction  according  to  the  magnitude  is  shown  in 
Tables  18,  19. 

Table  18.  The  Distribution  of  the  Terrain  Correction  According  to  tho 
Magnitude 

(5'  Mean  Block  Values) 


Interval  (mgals) 

Numbers 

Percentage 

0. 

.  10. 

9320682 

99.887  * 

10. 

.  20. 

5473 

0.059 

20. 

.  30. 

4652 

0.050 

30. 

.  40. 

321 

0.003 

40. 

.  50. 

30 

0.000 

50. 

.  60. 

9 

0.000 

60. 

.  70. 

3 

0.000 

70. 

.  80, 

6 

0.000 

80. 

.  90. 

9 

0.000 

90. 

.  100. 

9 

0.000 

100. 

.  110. 

2 

0.000 

110. 

.  120. 

2 

0.000 

120. 

.  130. 

0 

0.000 

Table  19.  The  Distribution  of  the  Terrain  Correction  According  to  the 
Magnitude 

(30'  Mean  Block  Values) 


Interval  (mgals) 

Numbers 

Percentage 

.  10. 

259140 

99.977  X 

.  20. 

5473 

0.022 

.  30. 

4652 

0.001 

The  statistics  of  the  gi  term  and  the  terrain  correction  are  exhibited  in 
Tables  20  and  21. 

Table  20.  Statistics  of  the  gi  Term  in  5’  and  30’  Mean  Block  Values 
Unit:  mgal 


Block  Size 

Mean  Value 

Standard  Dev. 

Min.  Value 

5’ 

0.27 

‘2.56 

442.14 

-78.88 

30’ 

0.27 

‘1.54 

45.08 

-10.47 

r 

0.27 

‘1.24 

25.52 

-5.16 

Table  21.  Statistics  of  the  Terrain  Correction  in  5*  and  30*  Mean  Block 
Values;  Unit:  mgal 


1 

Mean  Value 

Standard  Dev. 

Max.  Value 

Min.  Value 

5’ 

0.23 

‘1.01 

0.0 

30’ 

0.23 

‘0.-82 

0.0 

r 

0.23 

‘0174 

17.77 

0.0 

The  locations  of  the  absolute  values  of  the  gi  term  greater  than  5  mgals 
are  plotted  in  Figure  6.  The  locations  of  the  terrain  correction  greater  than  5 
mgals  are  shown  in  Figure  7. 
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LONGITUDE 
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Figure  6.  Location  of  g,  Values  that  Exceed  S  iiigul  in  30' x  30'  Blocks 


LONGITUDE 


Figure  7.  Location  of  Terrain  Correction  Valiiea  that  Exceed  5  tngul  in  30  x30 
Blocks 


1.  Compariaon  of  the  Degree  Variances  of  the  Gravity  Anomaly,  the  g.  Term 
and  the  Terrain  Correction 


The  corrections  of  the  disturbing  potential  caused  by  the  gi  term  and  the 
terrain  correction  can  be  expanded  in  spherical  harmonics  as  follows: 


dTi  = 

r 

N 

1 

n-2 

n 

1 

m-o 

(f) 

(C‘„,  cos  mX  +  S^n,  sin  mX) 

• 

l*nin 

(sin 

♦) 

i  =  1,2 

(66) 

where: 


(5T< 

kM 

_  _  a 

^nm  »  ^nm 
^nm 

r,  ♦,  X 


correction  of  the  disturbing  potential  caused  by  the  gi 

term  and  TC  for  i  =  1,  2,  respectively 

geocentric  gravitational  constant 

equational  radius  of  the  reference  ellipsoid 

fully  normalized  potential  coefficients 

fully  normalized  Legendre  function  of  degree  n  and 

order  m 

geocentric  coordinates 


The  degree  variance  of  the  disturbing  potential  is  defined  as  (Rapp, 
1986,  p.lO): 


n 


=  E  (CS„  + 

m-o 


(67) 


For  the  degree  variance,  c„,  of  the  gravity  anomaly  we  have: 

c„  =  7^  (n-1)^  <r„  (68) 

where  y  -  kM/a*.  In  the  same  manner  the  degree  variances  of  the 
corrections  to  the  potential  are: 


<  =  I  ((C'„)*  ^  (S‘,„)»]  i  =  1.2  (69) 

m—  o 

and  the  anomaly  degree  variances  of  the  g,  term  and  TC  are 

c^  =  7*  (n-1)*  ”  1.2  (70) 

In  order  to  get  an  impression  about  the  corrections  of  the  g,  term 
and  the  terrain  correction  on  the  gravity  anomaly  the  gravity  model 
0SUU6E  has  been  used.  The  percentage  of  the  root  mean  square  (RMS)  of 
the  degree  variances  is  defined  as  following 

Percentage  =  •  IQOX  i  =  1,2  (71) 

<fn 

The  gi  term  and  TC  in  5’  mean  block  values  were  used  to  form  I'x  1* 
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block  mean  values.  The  program  F419BV3  from  Rapp’s  program  library 
has  been  used  to  determine  the  coefficients  of  the  harmonics  spherical  of 
the  gi  term  and  the  TC  up  to  the  180*^^  degree. 

The  results  are  drawn  in  Figure  8.  The  RMS  values  of  the  degree  variances 
of  the  g,  term  and  TC  are  almost  2  percent  of  the  RMS  value  of  the  degree 
variance  of  0SU86E.  This  is  smaller  than  expected. 


o 

o 


DEGREES 


Figure  8.  Ratio  of  the  RMS  Values  of  the  Degree  Variances  of  the  gi  Term 
and  TC  to  the  Gravity  Model  0SU86E 

The  influence  of  the  topography  on  the  gravity  anomaly  and  its  spherical 
harmonic  expansion  was  studied  by  Rapp  (1977).  Different  models  of  the 
terrain  correction  were  used  and  the  correction  terms  found  were  on  the 
order  of  10  to  30%  at  the  lower  degree  of  the  spherical  harmonics  of  the 
gravity  anomaly.  The  reason  for  the  difference  should  be  due  to  adopted 
terrain  correction  models. 

It  is  very  interesting  to  find  that  the  coefficients  of  the  lower  degree  of 
the  g,  term  and  the  terrain  correction  are  almost  the  same.  Figure  9  shows 
the  ratio  of  the  RMS  of  variances  of  the  gi  term  to  TC. 
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o 

o 


DEGREES 


Figure  9.  Ratio  of  the  RMS  Value  of  the  Degree  Variance  of  the  gi  Term  to 
the  RMS  Value  of  the  Degree  Variance  of  TC 


Theoretically,  from  (12)  we  can  see  that  the  g,  term  equals  the 
terrain  correction  plus  the  integral 

-  -  KpR*  //  d<r  (72) 

The  computation  shows  this  integral  has  no  significant  contribution 
on  the  lower  degree  of  the  coefficients  of  the  spherical  harmonics  of 
the  gi  term. 

8.  The  Influence  of  the  g,  Term  and  TC  on  the  Geoid  Undulation  and 
the  Deflections  of  the  Vertical 

Now  we  want  to  know  the  influences  of  the  g]  term  and  TC  .on  the 
geoid  undulation  and  the  deflections  of  the  vertical. 

The  corrections  of  the  geoid  undulation  due  to  g,  term  and  TC  are 

«N*  =  —  i  =  1,2  .  (73) 

r 
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Using  (66^  and  (73)  we  have  the  mean  square  of  the  corrections  of 
the  geoid  undulation  due  to  the  g,  term  and  TC 

=  M  {(<5N<)»}  =  a*  !  <T«  i  =  1,2  ,  (74 

n=o 

where  M  {  }  denotes  the  average  over  the  whole  earth.  Its  definition 

is  given  in  (Heiskanen  and  Moritz,  1967,  p.252).  The  RMS  of  the  corrections 
of  the  geoid  undulation  is 


J  (6N‘)* 


i  =  1,2 


(75) 


The  vedues  of  equation  (75)  have  been  eveduated  for  the  gi  expansion  and 
the  terrain  correction  expansion.  These  values,  for  selected  degrees,  are 
given  in  Table  22. 

Table  22.  RMS  Value  of  the  Correction  Terms  to  the  Geoid  Undulations; 
Unit:  meter 


Degree 

gi  term 

Terrain  Correction 

2 

0.54 

0.54 

4 

0.13 

0.13 

5 

0.07 

0.07 

10 

0.04 

0.04 

15 

0.02 

0.02 

20 

0.01 

0.01 

30 

0.01 

0.01 

36 

0.01 

0.01 

50 

0.00 

0.00 

100 

0.00 

0.00 

150 

0.00 

0.00 

100 

0.00 

0.00 

We  have  two  conclusions  from  the  Table  22;  the  contributions  of  the  gi 
term  and  the  terrain  correction  to  the  geoid  undulation  are  the  same;  only  the 
lower  degrees  of  the  spherical  harmonic  expansions  of  the  correction  terms 
have  significant  effect  on  the  geoid  undulation.  The  effect  of  the  high 
degrees  (i>36)  harmonics  of  the  correction  terms  are  smaller  than  1 
cm. 


The  corrections  of  the  deflections  of  the  vertical  due  to  the  gi  term  and 
TC  are 

go*  =  -  grad  gN*  (76) 


where  0  is  the  vector  of  the  deflections  of  the  vertical;  it  has  two 
components ,  (  ,  t). 

Using  the  orthogonal  property  of  the  spherical  harmonics  (Heiskanen 
and  Moritz,  1967,  p.262)  we  have 
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(60‘)^  =  M  {(60*)*} 


(n+1)  <r^ 


i  =  1.2 


(77) 


with 

(69*)^  =  (6(*)^  +  (fit,*)* 

The  RMS  value  of  the  correction  to  the  deflections  of  the  vertical  is 

vWtJ  =  /J.  <  i  =  ‘-2  (7“> 

Table  23  shows  the  RMS  value  of  the  corrections  of  the  geoid  undulation 
and  the  deflections  of  the  vertical.  They  are  in  the  order  1  meter  and  0.1". 
Here  we  must  recognize  that  the  corrections  are  mean  values  the  whole  earth 
which  are  70%  covered  by  water  on  which  the  g,  term  and  TO  are 
zero. 

Table  23 .  RMS  Value  of  the  Correction  of  the  Geoid  Undulation  and  the 
Deflections  of  the  Vertical  Due  to  the  Term  and  TC 
Unit:  6N  in  iueter;  69  in  second 


correction 

RMS  of  6N 

RMS  of  66 

Si 

0.7071 

0.1129 

TC 

0.7065 

0.0880 

9.  Conclusion 


In  order  to  determine  the  coefficients  of  the  spherical  harmonic  expansion 
of  the  disturbing  potential  of  the  earth  in  a  more  precise  manner  the  gravity 

anomalies  have  to  be  analytically  downward  continued  to  a  simple  surface  - 

the  chosen  ellipsoid.  The  basic  formulas  for  the  downward  continuation  are 
the  Poisson's  integral  and  the  gradient  solution  (the  gj  term)  which  is  an 
approximation  of  the  Poisson’s  integral.  The  terrain  correction  has  been 
considered  as  an  approximate  method  to  eliminate  the  effect  of  the  topogiuphy 
above  the  ellipsoid. 

Because  we  have  no  dense  gravity  anomaly  data  on  a  global  basis,  we 
approximate  the  gravity  anomaly  by  using  the  assumption  (10).  The  elevation 
in  5’  mean  block  values  has  been  used  to  compute  the  g,  term  and  the  terrain 

correction  on  a  global  basis.  Before  we  started  the  work,  some  preliminary 

tests  were  done.  The  digital  elevation  model  in  30'  '  interval  in  the  United 
States  has  been  used  for  the  tests.  The  comparisons  of  the  results  which  use 
30'  '  point  elevation  and  the  elevation  in  5'  mean  block  values  have  been 
made.  The  results  show  that  the  computation  using  the  elevation  in  5'x  5' 
mean  block  values  are  acceptable. 

After  the  computations  of  the  g,  term  and  TC  in  5’  mean  block  values  in  a 

global  basis,  gi  and  TC  are  expanded  in  the  spherical  harmonica.  The 

influences  of  the  corrections  (g^  and  TC)  on  the  disturbing  potential  of  the 
earth  have  been  studied.  The  corrections  have  the  order  of  about  2%  of  the 
RMS  value  of  degree  variances  of  the  disturbing  potential.  It  is  much  smaller 
than  expected,  especially  at  the  lower  degree  of  the  spherical  harmonic  of  the 
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disturbing  potential.  The  influences  of  the  gi  term  and  TC  on  the  geoid 
undulation  (when  estimated  from  expansions  of  surface  gravity  data  only)  have 
been  considered.  It  takes  the  order  of  1  meter  (RMS  of  correction  of  the 
geoid  undulation).  For  the  deflections  of  the  vertical  the  RMS  of  the 
corrections  of  the  gi  term  and  TC  is  on  the  order  0.1". 

There  are  two  main  error  sources  in  the  computations:  the  assumption  (10) 
and  the  use  of  the  5’  mean  block  values.  They  both  are  due  to  lack  of  the 
denser  gravity  anomalies  and  denser  elevation  data  on  a  global  basis.  If  the 
denser  data  sets  are  available,  then  the  downward  continuation  of  the  gravity 
anomaly  to  the  ellipsoid  can  be  done  more  precisely.  We  can  then  hope  to  get 
a  more  accurate  downward  continuation  solution  needed  for  the  spherical 
harmonic  expansion  of  the  disturbing  potential  of  the  earth. 
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